#ifndef PLANT1D_SPDE_INPUT_1DSIN_H
#define PLANT1D_SPDE_INPUT_1DSIN_H

#include "Signal.h"

/**************************************************************
class plant2D_SPDE_input_1Dsin solve the following 2D EW equation
\[
\frac{\partial h}{\partial t} = w(x,y)+c_2\left(\frac{\partial^2 h}{\partial x^2}+\frac{\partial^2 h}{\partial y^2}\right)+\xi(x,y,t)
\]
where $x\in[0,XYmax],y\in[0,XYmax]$. $\xi(x,y,t)$ is a Gaussian
white noise with zero mean and the following covariance
\[
\left<\xi(x,y,t)\xi(x',y',t')\right>=\sigma^2\delta(x-x')\delta(y-y')\delta(t-t').
\]
$w$ is a patterned deposition rate profile with the form
\[
w(x,y) = W+A\sin(\frac{2*k*\pi}{L}).
\]
Boundary condition:
h(0,y,t) = h(L,y,t),    h(x,0,t) = h(x,L,t)
\frac{\partial h}{\partial x}(0,y,t) = \frac{\partial h}{\partial x}(L,y,t)
\frac{\partial h}{\partial y}(x,0,t) = \frac{\partial h}{\partial y}(x,L,t)

Initial condition:
h(x,y,0) = h_0(x,y) = 0

In the code:
XYmax = L
Freq = k

class plant2D_SPDE_input_1Dsin use model decomposition method to
solvet the EW equation. Detailed can be found elsewhere.
**************************************************************/

class plant2D_SPDE_input_1Dsin:public block{
  // Inputs =====================================
  double T;        // unit: K
  double W;        // unit: ?
  double A;        // unit: ?
  double Freq;     // unit: number of periods in the domain

  // Model parameters ===========================
  int LatticeSize; // number of discrete points in the h profile
  double Xmax;    // unit: nm

  // Numerical solution parameters =============
  int mode;        // number of mode used when solving EW equation 
  // For analytical expression
  double sqrt_dt;
  // For explicit/implicit method
  // block::dt is the time interval for output, dt_con is the dt used in numerical integration. 
  double dt_con;
  double sqrt_dt_con;
  double dt_conservative_factor;

  // States =====================================
  double *z1;
  double *z2;
  double *z1_next;
  double *z2_next;
  
  // Model parameters ===========================
  double *lambda;
  double *Km2[2];
  
  // Random number generator=====================
  const gsl_rng_type * noise_type;
  gsl_rng * noise;

  double eq_c2(double W);
  double eq_sigma2(double W);

  inline int index1D(int m){return m;};
  // Eigenfunctions =============================
  inline double phi1m(int m,int ix){
    return 2.0*sin(2*m*pi*ix/LatticeSize)/Xmax;
  };
  
  inline double phi2m(int m,int ix){
    double coeff;
    if(m == 0 || m*2==LatticeSize){
      coeff = 1.0/Xmax;
    }
    else if(m!=0 && 2*m!=LatticeSize){
      coeff = 2.0/Xmax;
    }
    else{
      coeff = SQRT_TWO/Xmax;
    }
    return coeff*cos(2*m*pi*ix/LatticeSize);
  };

  

  
 public:
  plant2D_SPDE_input_1Dsin(int M,int L,double idt,double iXYmax);
  ~plant2D_SPDE_input_1Dsin();
  virtual void update(double sysTime);
  virtual void reset();
  void output(string filename);
  double getM2(); // Return m^2, unit = dimensionless
  double getR2(); // Return r^2, unit = nm^2
  double getH();  // Return hmean, unit = nm
  double getCorLength(); // Return correlation length, unit = nm
  void   getSnapshot(string FileName);
};
#endif
